#include "mcmcModelSupervisor.h"

MCMC_Model* CreateModel(const int type)
{
	MCMC_Model *pmodel=NULL;

	switch(type)
	{
	case MODEL_LUVCOLOR_GAUSS:
		pmodel = new MCMC_LuvColor_Gaussian();		
		break;
	case MODEL_LUVCOLOR_MIXTUREGAUSS:
		pmodel = new MCMC_LuvColor_MixtureGaussians();		
		break;
	case MODEL_2DCOLOR_GAUSS:
		pmodel = NULL;
		assert(false);
		break;
	case MODEL_3DCOLOR_GAUSS:
		pmodel = NULL;
		assert(false);
		break;
	case MODEL_3DCOLOR_MIXGAUSS:
		pmodel = NULL;
		assert(false);
		break;
	case MODEL_TEXTURE:
		pmodel = NULL;
		assert(false);
		break;
	case MODEL_INTENSITY_GAUSS:
		pmodel = new MCMC_Intensity_Gaussian();
		break;
	case MODEL_INTENSITY_MIXGAUSS:
		pmodel = NULL;
		assert(false);
		break;
	case MODEL_INTENSITY_TEXTURE:
		pmodel = NULL;
		assert(false);
		break;
	case MODEL_INTENSITY_DELTA:
		pmodel = new MCMC_Intensity_Delta();
		break;
	case MODEL_INTENSITY_SURFACE:
		pmodel = new MCMC_Intensity_Surface();
		break;
	case MODEL_LUVCOLOR_SURFACE:
		pmodel = new MCMC_LuvColor_Surface();		
		break;
	}
	return pmodel;

	//MCMC_Model *pmodel=NULL;

	//switch(type)
	//{
	//case MODEL_LUVCOLOR_GAUSS:
	//	pmodel = new MCMC_LuvColor_Gaussian();
	//	break;
	//case MODEL_LUVCOLOR_MIXTUREGAUSS:
	//	pmodel = new MCMC_LuvColor_MixtureGaussians();
	//	break;
	//case MODEL_2DCOLOR_GAUSS:
	//	pmodel = new MCMC_2DColor_Gaussian();
	//	break;
	//case MODEL_3DCOLOR_GAUSS:
	//	pmodel = new MCMC_3DColor_Gaussian();
	//	break;
	//case MODEL_3DCOLOR_MIXGAUSS:
	//	pmodel = new MCMC_3DColor_MixtureGaussians();
	//	break;
	//case MODEL_TEXTURE:
	//	pmodel = new MCMC_Model_Texture();
	//	break;
	//case MODEL_INTENSITY_GAUSS:
	//	pmodel = new MCMC_Intensity_Gaussian();
	//	break;
	//case MODEL_INTENSITY_MIXGAUSS:
	//	pmodel = new MCMC_Intensity_MixtureGaussians();
	//	break;
	//case MODEL_INTENSITY_TEXTURE:
	//	pmodel = new MCMC_Intensity_Texture();
	//	break;
	//case MODEL_INTENSITY_DELTA:
	//	pmodel = new MCMC_Intensity_Delta();
	//	break;
	//case MODEL_INTENSITY_SURFACE:
	//	pmodel = new MCMC_Intensity_Surface();
	//	break;
	//case MODEL_LUVCOLOR_SURFACE:
	//	pmodel = new MCMC_LuvColor_Surface();
	//	break;
	//}
	//return pmodel;
}


bool CopyModel(MCMC_Model *pmodel1, MCMC_Model *pmodel2)
{
	bool bResult=false;

	if (pmodel1!=NULL && pmodel2!=NULL && pmodel1->Type()==pmodel2->Type())
	{
		bResult = true;

		switch(pmodel1->Type())
		{
		case MODEL_LUVCOLOR_GAUSS:
			*((MCMC_LuvColor_Gaussian *)pmodel1)=*((MCMC_LuvColor_Gaussian *)pmodel2);
			break;
		case MODEL_LUVCOLOR_MIXTUREGAUSS:
			*((MCMC_LuvColor_MixtureGaussians *)pmodel1)=*((MCMC_LuvColor_MixtureGaussians *)pmodel2);
			break;
		case MODEL_TEXTURE:
			assert(false);
			break;
		case MODEL_INTENSITY_GAUSS:
			*((MCMC_Intensity_Gaussian *)pmodel1)=*((MCMC_Intensity_Gaussian *)pmodel2);
			break;
		case MODEL_INTENSITY_TEXTURE:
			assert(false);
			break;
		case MODEL_INTENSITY_DELTA:
			*((MCMC_Intensity_Delta *)pmodel1)=*((MCMC_Intensity_Delta *)pmodel2);
			break;
		case MODEL_INTENSITY_SURFACE:
			*((MCMC_Intensity_Surface *)pmodel1)=*((MCMC_Intensity_Surface *)pmodel2);
			break;
		case MODEL_LUVCOLOR_SURFACE:
			*((MCMC_LuvColor_Surface *)pmodel1)=*((MCMC_LuvColor_Surface *)pmodel2);
			break;
		default:
			assert(false);
		}
	}
	return bResult;

	/*bool bResult=false;

	if (pmodel1!=NULL && pmodel2!=NULL && pmodel1->Type()==pmodel2->Type())
	{
		bResult = true;

		switch(pmodel1->Type())
		{
		case MODEL_LUVCOLOR_GAUSS:
			*((MCMC_LuvColor_Gaussian *)pmodel1)=*((MCMC_LuvColor_Gaussian *)pmodel2);
			break;
		case MODEL_LUVCOLOR_MIXTUREGAUSS:
			*((MCMC_LuvColor_MixtureGaussians *)pmodel1)=*((MCMC_LuvColor_MixtureGaussians *)pmodel2);
			break;
		case MODEL_2DCOLOR_GAUSS:
			*((MCMC_2DColor_Gaussian *)pmodel1)=*((MCMC_2DColor_Gaussian *)pmodel2);
			break;
		case MODEL_3DCOLOR_GAUSS:
			*((MCMC_3DColor_Gaussian *)pmodel1)=*((MCMC_3DColor_Gaussian *)pmodel2);
			break;
		case MODEL_3DCOLOR_MIXGAUSS:
			*((MCMC_3DColor_MixtureGaussians *)pmodel1)=*((MCMC_3DColor_MixtureGaussians *)pmodel2);
			break;
		case MODEL_TEXTURE:
			*((MCMC_Model_Texture *)pmodel1)=*((MCMC_Model_Texture *)pmodel2);
			break;
		case MODEL_INTENSITY_GAUSS:
			*((MCMC_Intensity_Gaussian *)pmodel1)=*((MCMC_Intensity_Gaussian *)pmodel2);
			break;
		case MODEL_INTENSITY_MIXGAUSS:
			*((MCMC_Intensity_MixtureGaussians *)pmodel1)=*((MCMC_Intensity_MixtureGaussians *)pmodel2);
			break;
		case MODEL_INTENSITY_TEXTURE:
			*((MCMC_Intensity_Texture *)pmodel1)=*((MCMC_Intensity_Texture *)pmodel2);
			break;
		case MODEL_INTENSITY_DELTA:
			*((MCMC_Intensity_Delta *)pmodel1)=*((MCMC_Intensity_Delta *)pmodel2);
			break;
		case MODEL_INTENSITY_SURFACE:
			*((MCMC_Intensity_Surface *)pmodel1)=*((MCMC_Intensity_Surface *)pmodel2);
			break;
		case MODEL_LUVCOLOR_SURFACE:
			*((MCMC_LuvColor_Surface *)pmodel1)=*((MCMC_LuvColor_Surface *)pmodel2);
			break;
		}
	}
	return bResult;*/
}

//////////////////////////////////////////////////////
//////////////////////////////////////////////////////
// construction & destruction
MCMC_ModelSupervisor::MCMC_ModelSupervisor()
{
	pmap				= NULL;
	pparticle_gauss		= NULL;
	prespFilter			= NULL;
}

MCMC_ModelSupervisor::~MCMC_ModelSupervisor()
{
	pmap				= NULL;
	pparticle_gauss		= NULL;
	prespFilter			= NULL;
}

MCMC_Particle* MCMC_ModelSupervisor::Particle(const int model_type)
{
	if (model_type==MODEL_TEXTURE || model_type==MODEL_INTENSITY_TEXTURE)
	{
		assert(false);
		return NULL;
	}
	else if (model_type==MODEL_INTENSITY_DELTA || model_type==MODEL_INTENSITY_SURFACE ||
		     model_type==MODEL_LUVCOLOR_SURFACE)
		return NULL;
	else
		return pparticle_gauss;
}

void MCMC_ModelSupervisor::Set(MCMC_Map_2Layers &map, MCMC_Particle &particle1,MCMC_FilterResp &respFilter, const MCMC_SegSetting& setting)
{
	pmap			= &map;
	pparticle_gauss	= &particle1;
	prespFilter		= &respFilter;
	bcolor			= setting.bcolor;

	double qs[5];

	// set up proposal probabilities for different families of models
	if (setting.btexture)
	{
		qs[0]	= 0.45;		// Gaussian
		qs[1]	= 0.00;		// Mixture of Gaussians
		qs[2]	= 0.20;		// Delta
		qs[3]	= 0.20;		// Texture
		qs[4]	= 0.15;		// Surface
	}
	else
	{
		if (!IsColorImage())
		{
			// gray scale image
			qs[0]	= 0.55;		// Gaussian
			qs[1]	= 0.00;		// Mixture of Gaussians
			qs[2]	= 0.25;		// Delta
			qs[3]	= 0.00;		// Texture
			qs[4]	= 0.20;		// Surface
		}
		else
		{
			// color image
			qs[0]	= 0.65;		// Gaussian
			qs[1]	= 0.00;		// Mixture of Gaussians
			qs[2]	= 0.00;		// Delta
			qs[3]	= 0.00;		// Texture
			qs[4]	= 0.15;		// Surface
		}
	}

	double sum=0.0;
	for (int i=0; i<5; i++)		sum += qs[i];
	for (int i=0; i<5; i++)			qs[i] = qs[i]/sum;

	q_model_gauss		= qs[0];
	q_model_mixgauss	= qs[1];
	q_model_delta		= qs[2];
	q_model_texture		= qs[3];
	q_model_surface		= qs[4];

	// set up bottom-up cues
	bottomup_cues.Set(respFilter);
}

void MCMC_ModelSupervisor::InitRegionModels(void)
{
	MCMC_RegionModel	*pmodel;

	pmodel = NULL;
	// setup the scale factor
	scale_factor = LayerRegion()->scale_factor;

	// initialization
	LinkList<int>	list_topnode_regions,list_aregions;
	bool			bgo;
	int				indexid;
	MCMC_RegionElement *pregion;

	LayerRegion()->GetAllTopNodeRegions(list_topnode_regions);
	bgo = list_topnode_regions.GetFirst(indexid);

	
	while (bgo)
	{
		pregion = LayerRegion()->GetRegionFromIndexId(indexid);
		
		assert(pregion!=NULL);		

		if (pmodel != NULL)
			delete pmodel;

		pregion->GetAllAtomicRegions(list_aregions);        
		pmodel = ProposeModel(list_aregions);		
		pregion->SetStat(pmodel);
	
		//// shape prior
		pregion->ComputeEnergyShape();
		
		bgo = list_topnode_regions.GetNext(indexid);
		
	}
	
	if (pmodel!=NULL)
		delete pmodel;
	
}

double MCMC_ModelSupervisor::GetPdf(MCMC_Model *pmodel)
{
	double	P,fraction;

	if (pmodel == NULL)
		return 0.0;

	switch(pmodel->Type())
	{
	case MODEL_INTENSITY_GAUSS:
	case MODEL_3DCOLOR_GAUSS:
	case MODEL_LUVCOLOR_GAUSS:
		fraction=((double)pmodel->vc_vote.GetMax())/pmodel->vc_vote.Sum();
		if (fraction>0.15 && fraction<0.85)
			P = q_model_gauss;
		else
			P = q_model_gauss+q_model_mixgauss;
		break;
	case MODEL_INTENSITY_MIXGAUSS:
	case MODEL_3DCOLOR_MIXGAUSS:
	case MODEL_LUVCOLOR_MIXTUREGAUSS:
		fraction=((double)pmodel->vc_vote.GetMax())/pmodel->vc_vote.Sum();
		if (fraction>0.15 && fraction<0.85)
			P = q_model_mixgauss;
		else
			P = 0.0;
		break;
	case MODEL_INTENSITY_DELTA:
		P = q_model_delta;
		break;
	case MODEL_TEXTURE:
	case MODEL_INTENSITY_TEXTURE:
		if (pmodel->nsize<1700)
			P = 0.0;
		else
			P = q_model_texture;
		break;
	case MODEL_INTENSITY_SURFACE:
	case MODEL_LUVCOLOR_SURFACE:
		if (pmodel->nsize<400)
			P = 0.0;
		else
			P = q_model_surface;
		break;
	default : P = 1.0;
	}

	switch(pmodel->Type())
	{
	case MODEL_INTENSITY_GAUSS:
	case MODEL_INTENSITY_MIXGAUSS:
	case MODEL_INTENSITY_DELTA:
	case MODEL_INTENSITY_TEXTURE:
	case MODEL_INTENSITY_SURFACE:
		if (IsColorImage())
			P = 0.0;
		break;

	case MODEL_3DCOLOR_GAUSS:
	case MODEL_LUVCOLOR_GAUSS:
	case MODEL_3DCOLOR_MIXGAUSS:
	case MODEL_LUVCOLOR_MIXTUREGAUSS:
	case MODEL_LUVCOLOR_SURFACE:
		if (!IsColorImage())
			P = 0.0;
		break;
	}

	return pmodel->GetPdf()*P;
}

MCMC_Model*	MCMC_ModelSupervisor::SelectModelType(mcmcMatrix<int> &mx_label, const int nlabel, irmm_Rect &bound)
{
	MCMC_Model	*pmodel;
	int			imodel;
	RF			rf;

	// initialization
	rf.Set(1,NUM_MODEL_TYPE,1);
	for (imodel=1; imodel<=NUM_MODEL_TYPE; imodel++)
		rf.SetOneBin(imodel,0.0);

	if (!IsColorImage())
	{
		// gray scale image
		rf.SetOneBin(MODEL_INTENSITY_GAUSS, q_model_gauss);
		rf.SetOneBin(MODEL_INTENSITY_MIXGAUSS, q_model_mixgauss);
		rf.SetOneBin(MODEL_INTENSITY_DELTA, q_model_delta);
		rf.SetOneBin(MODEL_INTENSITY_TEXTURE, q_model_texture);
		rf.SetOneBin(MODEL_INTENSITY_SURFACE, q_model_surface);
	}
	else
	{
		// color image
#ifdef LUV_COLOR_MODEL
		rf.SetOneBin(MODEL_LUVCOLOR_GAUSS,q_model_gauss);
		rf.SetOneBin(MODEL_LUVCOLOR_MIXTUREGAUSS,q_model_mixgauss);
		rf.SetOneBin(MODEL_LUVCOLOR_SURFACE,q_model_surface);
#else
		rf.SetOneBin(MODEL_3DCOLOR_GAUSS, q_model_gauss);
		rf.SetOneBin(MODEL_3DCOLOR_MIXGAUSS, q_model_mixgauss);
#endif
	}

	// compute the number of pixels first
	int	count=0;
	double	grad=0.0;

	for (int j=bound.top; j<=bound.bottom; j++)
		for (int i=bound.left; i<=bound.right; i++)
		{
			if (mx_label(j,i)>0 && mx_label(j,i)==nlabel)
			{
				count++;
				grad += prespFilter->ras_edge_g(j,i);
			}
		}
	if (count<1700)
	{
		rf.SetOneBin(MODEL_TEXTURE,0.0);
		rf.SetOneBin(MODEL_INTENSITY_TEXTURE,0.0);
	}
	if (grad/((double)count)>2.0 || count<400)
	{
		rf.SetOneBin(MODEL_INTENSITY_SURFACE,0.0);
		rf.SetOneBin(MODEL_LUVCOLOR_SURFACE,0.0);
	}

	imodel = (int)rf.SamplingNonSetCdf();

	if (imodel!=MODEL_TEXTURE && imodel!=MODEL_INTENSITY_TEXTURE && imodel!=MODEL_INTENSITY_DELTA &&
		imodel!=MODEL_INTENSITY_SURFACE && imodel!=MODEL_LUVCOLOR_SURFACE)
	{
		mcmcVector<float> vc_vote;
		vc_vote.SetRowVector(pparticle_gauss->Num());
		vc_vote.InitValue(0);

		for (int j=bound.top; j<=bound.bottom; j++)
			for (int i=bound.left; i<=bound.right; i++)
			{
				if (mx_label(j,i)>0 && mx_label(j,i)==nlabel)
					pparticle_gauss->AddClusterVote(vc_vote,j,i);
			}
		double	fraction=((double)vc_vote.GetMax())/vc_vote.Sum();
		if (fraction<=0.15 || fraction>=0.85)
		{
			if (IsColorImage())
			{
#ifdef LUV_COLOR_MODEL
				imodel = MODEL_LUVCOLOR_GAUSS;
#else
				imodel = MODEL_3DCOLOR_GAUSS;
#endif
			}
			else
				imodel = MODEL_INTENSITY_GAUSS;
		}
	}
	pmodel = CreateModel(imodel);
	assert(pmodel != NULL);

	pmodel->SetImage(*prespFilter);
	((MCMC_RegionModel*) pmodel)->Set(mx_label,bound,*Particle(pmodel->Type()), NULL);

	return pmodel;
}

void MCMC_ModelSupervisor::ChooseModel(mcmcMatrix<int> &mx_label, const int num_label, irmm_Rect &bound, MCMC_Model* pmodels[])
{
	int				nindex;
	MCMC_Model		*pmodel;

	// select a cluster
	for (nindex=0; nindex<num_label; nindex++)
	{
		pmodel = SelectModelType(mx_label,nindex+1,bound);

		pmodel->Mean(nindex+1);

		pmodel->DiffuseComputedMean(nindex+1);
		// compute variance
		pmodel->Variance(nindex+1);
		// compute statistics
		pmodel->LogLikelihood(nindex+1);

		pmodels[nindex] = pmodel;
	}
}

void MCMC_ModelSupervisor::FitAllRegionParms(void)
{
	LinkList<int>		list_topnode_regions;
	bool				bgo;
	int					indexid;
	MCMC_RegionElement	*pregion;
	LinkList<PMCMC_PixelSet> list_domain;
	MCMC_RegionModel	*pmodel;

	LayerRegion()->GetAllTopNodeRegions(list_topnode_regions);
	bgo = list_topnode_regions.GetFirst(indexid);
	while (bgo)
	{
		pregion = LayerRegion()->GetRegionFromIndexId(indexid);

		pregion->GetRegionDomain(list_domain);
		
		assert(pregion->Model()!= NULL);
		
		pmodel = (MCMC_RegionModel*)pregion->Model();

		pmodel->SetImage(*prespFilter);
		
		pmodel->Set(list_domain, *Particle(pmodel->Type()), NULL);

		pmodel->Mean(0);

		pmodel->DiffuseComputedMean(0);
		// compute variance
		pmodel->Variance(0);
		// compute statistics
		pmodel->LogLikelihood(0);
		
		bgo = list_topnode_regions.GetNext(indexid);
	}
}



void MCMC_ModelSupervisor::RecomputePosteriorRegions(LinkList<int> &list_region_ids)
{
	int			id_region;
	MCMC_AtomicRegion *pregion;
	irmm_Rect	bound;
	bool		bgo;

	bgo = list_region_ids.GetFirst(id_region);
	while (bgo)
	{
		pregion = LayerRegion()->GetAtomicRegion(id_region);
		assert(pregion != NULL);

		bound = pregion->GetBound();

		pregion->Model()->SetMatrixLabel(pmap->ras_region_label);
		pregion->Model()->SetBound(bound);

		pregion->Model()->LogLikelihood(id_region);

		bgo = list_region_ids.GetNext(id_region);
	}
}

void MCMC_ModelSupervisor::RecomputePosteriorRegions(LinkList<int> &list_region_ids, MCMC_PixelSet &list_pixels_old,
											  MCMC_PixelSet &list_pixels_new)
{
	if (list_pixels_old.Num()<=0 && list_pixels_new.Num()<=0)
		return;		// no change

	bool		bgo,bfirst=false;
	int			id_region;
	MCMC_AtomicRegion		*pregion;

	bgo = list_region_ids.GetFirst(id_region);
	while (bgo)
	{
		pregion = LayerRegion()->GetAtomicRegion(id_region);
		assert(pregion != NULL);

		if (pregion->Valid() && pregion->Model()!=NULL)
		{
			pregion->Model()->SetMatrixLabel(pmap->ras_region_label);
			pregion->Model()->UpdateLogLikelihood(id_region, list_pixels_old, list_pixels_new);
		}

		bgo = list_region_ids.GetNext(id_region);
	}
}


// compute q(l|Tst(I))
/*
	TOC_Intensity=0,
	TOC_Smooth=1,
	TOC_Texture=2,
	TOC_Face=3,
	TOC_Text=4
*/
int MCMC_ModelSupervisor::PrepareModelType(RF& rf, LinkList<int> &list_aregions, LinkList<int> *plist_types_to_consider)
{

	//return num_non_zero;
	rf.Set(0, NUM_MODEL_TYPE, 1);
	rf.InitValue(0.0);

	// compute each individual cues

	bool	bgo1,bgo2; //,p_sign[NUM_MODEL_TYPE+1];
	int		id_aregion,j,i;
	int	num;
	double	dmean_overall;

	MCMC_AtomicRegion		*paregion;

	for (i=0; i<=NUM_MODEL_TYPE; i++)
		p_sign[i] = false;

	if (plist_types_to_consider == NULL)
	{
		/*for (i=0; i<=NUM_MODEL_TYPE; i++)
			p_sign[i] = false;*/

		if ( !IsColorImage() )
		{
			p_sign[TOC_Intensity] = true;
			p_sign[TOC_Smooth] = true;
			p_sign[TOC_Texture] = true;			
		}
		else
		{
#ifdef LUV_COLOR_MODEL
			p_sign[TOC_LUVColor] = true;
			p_sign[TOC_LUVColorMix] = true;
			p_sign[TOC_LUVSurface] = true;
#else
			p_sign[TOC_3DColor] = true;
			p_sign[TOC_3DColor_Mix] = true;
#endif
			
		}		
		p_sign[TOC_Text] = false;
		p_sign[TOC_Face] = false;
	}
	else
	{
		for (i=0; i<=NUM_MODEL_TYPE; i++)
			if (plist_types_to_consider->IsIn(i))
				p_sign[i] = true;
			else
				p_sign[i] = false;
	}

	// check for gaussian model
	if ( p_sign[TOC_Intensity] || p_sign[TOC_3DColor] || p_sign[TOC_LUVColor] )
	{
		// compute mean first
		dmean_overall = 0.0;
		num = 0;

		// retrieve pixels
		bgo1 = list_aregions.GetFirst(id_aregion);
		while (bgo1)
		{
			paregion = LayerRegion()->GetAtomicRegion(id_aregion);
			if (paregion != NULL)
			{
				bgo2 = paregion->GetFirstPixel(j, i);
				while (bgo2)
				{
					dmean_overall += prespFilter->ImageIntensity(j, i);
					num++;

					bgo2 = paregion->GetNextPixel(j, i);
				}
			}
			bgo1 = list_aregions.GetNext(id_aregion);
		}

		if (num > 0)
			dmean_overall /= (double)num;
		else
			dmean_overall = 0.0;
	}

	double dvar,dtemp,dgrad;
	MCMC_CueLocalGaussian	*pcue_gaussain;
	MCMC_CueEdge			*pcue_edge;

	num = 0;
	dvar = 0.0;
	dgrad = 0.0;
	// check for the model type
	// retrieve pixels
	bgo1 = list_aregions.GetFirst(id_aregion);
	while (bgo1)
	{
		paregion = LayerRegion()->GetAtomicRegion(id_aregion);
		if (paregion != NULL)
		{
			bgo2 = paregion->GetFirstPixel(j, i);
			while (bgo2)
			{
				if ( p_sign[TOC_Intensity] || p_sign[TOC_3DColor] || p_sign[TOC_3DColor_Mix] || p_sign[TOC_LUVColor] || p_sign[TOC_LUVColorMix] )
				{
					// gaussian model
					pcue_gaussain = (MCMC_CueLocalGaussian*) bottomup_cues.mx_cue_pixel(j, i).GetCue(0);
					dtemp = (pcue_gaussain->mean-dmean_overall);
					dvar += dtemp*dtemp;
				}

				if (p_sign[TOC_Smooth] || p_sign[TOC_Texture] || p_sign[TOC_LUVSurface])
				{
					pcue_edge = (MCMC_CueEdge*) bottomup_cues.mx_cue_pixel(j, i).GetCue(1);
					dgrad += pcue_edge->g;
				}

				num++;

				bgo2 = paregion->GetNextPixel(j, i);
			}
		}
		bgo1 = list_aregions.GetNext(id_aregion);
	}

	if (num <= 0)
		return 0;

	// propose label type
	double q_l;

	// gaussian
	if (p_sign[TOC_Intensity] || p_sign[TOC_3DColor] || p_sign[TOC_3DColor_Mix] || p_sign[TOC_LUVColor] || p_sign[TOC_LUVColorMix])
	{
		dvar = sqrt(dvar/((double)num))/40.0;
		q_l = 1.0-dvar;
		q_l = max(q_l, 0.3);
		q_l = min(q_l, 0.8);
		if (p_sign[TOC_Intensity]) rf.SetOneBin(TOC_Intensity, q_l);
		if (p_sign[TOC_3DColor]) rf.SetOneBin(TOC_3DColor, q_l);
		if (p_sign[TOC_3DColor_Mix]) rf.SetOneBin(TOC_3DColor_Mix, q_l);
		if (p_sign[TOC_LUVColor]) rf.SetOneBin(TOC_LUVColor, q_l);
		if (p_sign[TOC_LUVColorMix]) rf.SetOneBin(TOC_LUVColorMix, q_l);
	}

	// smooth
	if (p_sign[TOC_Smooth] || p_sign[TOC_LUVSurface] )
	{
		q_l = 1.0-dgrad/((double)num);
		if (num < 50)
			q_l = 0.0;
		else
		{
			q_l = max(q_l, 0.3);
			q_l = min(q_l, 0.8);
		}

		if (p_sign[TOC_Smooth]) rf.SetOneBin(TOC_Smooth, q_l);
		if (p_sign[TOC_LUVSurface]) rf.SetOneBin(TOC_LUVSurface, q_l);
	}

	// texture
	if (p_sign[TOC_Texture])
	{
		q_l = dgrad/((double)num);

		q_l = max(q_l, 0.2);
		q_l = min(q_l, 0.8);

		rf.SetOneBin(TOC_Texture, q_l);
	}

	int num_non_zero=0;
	rf.normalize();

	for (i=0; i<=NUM_MODEL_TYPE; i++)
	{
		if (rf.pdf(i) > 0.0)
		{
			num_non_zero++;
		}
	}

	return num_non_zero;
}

double MCMC_ModelSupervisor::ProposeModelType(int &imodel, LinkList<int> &list_aregions, LinkList<int> *plist_types_to_consider)
{
	RF rf;
	int num_model;

	num_model = PrepareModelType(rf, list_aregions, plist_types_to_consider);
	if (num_model <= 0)
	{
		imodel = -1;
		return 0.0;
	}
	else
	{
		int toc;

		toc = IntClose(rf.SamplingNonSetCdf());

		imodel = ToModelType(toc);
		return rf.pdf(toc);
	}
}

double MCMC_ModelSupervisor::QProposeModelType(int imodel, LinkList<int> &list_aregions, LinkList<int> *plist_types_to_consider)
{
	int toc;

	toc = ToCueType(imodel);
	RF rf;
	PrepareModelType(rf, list_aregions, plist_types_to_consider);
	return rf.pdf(toc);
}


int	MCMC_ModelSupervisor::ToModelType(int toc)
{
	switch (toc)
	{
	case TOC_Intensity: return MODEL_INTENSITY_GAUSS;
	case TOC_Smooth: return MODEL_INTENSITY_SURFACE;
	case TOC_Texture: return MODEL_INTENSITY_DELTA;
	case TOC_Face: return MODEL_INTENSITY_FACE;
	case TOC_Text: return MODEL_INTENSITY_TEXT;
	case TOC_2DColor: return MODEL_2DCOLOR_GAUSS;
	case TOC_3DColor: return MODEL_3DCOLOR_GAUSS;
	case TOC_3DColor_Mix: return MODEL_3DCOLOR_MIXGAUSS;
	case TOC_LUVColor: return MODEL_LUVCOLOR_GAUSS;
	case TOC_LUVColorMix: return  MODEL_LUVCOLOR_MIXTUREGAUSS;
	case TOC_LUVSurface: return MODEL_LUVCOLOR_SURFACE;
	}
	return -1;
}

int	MCMC_ModelSupervisor::ToCueType(int imodel)
{
	switch (imodel)
	{
	case MODEL_INTENSITY_GAUSS: return TOC_Intensity;
	case MODEL_INTENSITY_SURFACE: return TOC_Smooth;
	case MODEL_INTENSITY_DELTA: return TOC_Texture;
	case MODEL_INTENSITY_FACE: return TOC_Face;
	case MODEL_INTENSITY_TEXT: return TOC_Text;
	case MODEL_2DCOLOR_GAUSS: return TOC_2DColor;
	case MODEL_3DCOLOR_GAUSS: return TOC_3DColor;
	case MODEL_3DCOLOR_MIXGAUSS: return TOC_3DColor_Mix;
	case MODEL_LUVCOLOR_GAUSS: return TOC_LUVColor;
	case MODEL_LUVCOLOR_MIXTUREGAUSS: return TOC_LUVColorMix;
	case MODEL_LUVCOLOR_SURFACE: return TOC_LUVSurface;
	}
	return -1;
}

MCMC_RegionModel* MCMC_ModelSupervisor::ProposeModel(LinkList<int> &list_aregions, LinkList<int> *plist_types_to_consider)
{
	int imodel;

	// propose a model type first
	ProposeModelType(imodel, list_aregions, plist_types_to_consider);
	
	return CreateAndComputeModel(imodel, list_aregions);
}



MCMC_RegionModel* MCMC_ModelSupervisor::CreateAndComputeModel(int imodel, LinkList<int> &list_aregions)
{
	MCMC_RegionModel *pmodel;
	
	pmodel = (MCMC_RegionModel *)CreateModel(imodel);
	assert(pmodel != NULL);
	pmodel->SetImage(*prespFilter);

	// pixel domain
	LinkList<PMCMC_PixelSet> list_domain;
	MCMC_AtomicRegion *paregion;
	int					id_aregion;
	bool bgo;

	if (list_aregions.Num() == 1)
	{
		assert(list_aregions.GetFirst(id_aregion));
		paregion = LayerRegion()->GetAtomicRegion(id_aregion);
		assert(paregion != NULL);

		list_domain.Insert(&paregion->pixel_set);
	}
	else
	{
		bgo = list_aregions.GetFirst(id_aregion);
		while (bgo)
		{
			paregion = LayerRegion()->GetAtomicRegion(id_aregion);
			if (paregion != NULL)
				list_domain.Insert(&paregion->pixel_set);
			
			bgo = list_aregions.GetNext(id_aregion);
		}
	}

	pmodel->Set(list_domain, *Particle(pmodel->Type()), NULL);
	
	pmodel->Mean(0);

	pmodel->DiffuseComputedMean(0);
	// compute variance
	pmodel->Variance(0);
	// compute statistics
	pmodel->LogLikelihood(0);

	return pmodel;
}
